% Copyright (c) 2012, John Gebbie. All rights reserved.
%
% This library is free software; you can redistribute it and/or modify it
% under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation; either version 2.1 of the License, or (at
% your option) any later version.
%
% This library is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
% General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this library; if not, write to the Free Software Foundation,
% Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA


function [ spect_full ] = fftsynth(spect_partial,freq_low,freqs_full,dim)
% [ spect_full ] = fftsynth(spect_partial,freq_low,freqs_full,dim)
%
%   Perform Fourier synthesis. This function expands signal specified with
%   a subset of the positive frequencies to include the full set of
%   positive and negative frequencies.  Negative frequency bins are
%   populated with the conjugate of the positive frequencies. The output is
%   of this function can be passed to ifft to produce a real time-domain
%   signal.
%
% INPUTS
%       spect_partial   : input signals (selected contiguous bins)
%       freq_low        : lowest frequency in spect_partial. if specified
%                           as a vector, only the first element is used.
%       freqs_full      : full set of frequency bins to output (see
%                           fftinfo)
%       dim             : dimension along which signals are laid out. if
%                           signals are in columns, dim is 1. default is 1.
%
% OUTPUTS
%       spect_full     : signal including full set of frequency bins
%
% EXAMPLE
%
%     fs = 1000; % sampling frequency (Hz)
%     fftsize = 250; % number of samples
%     fr_band = [100 200];
%     tau = 0.1;
%     fr_full = fftinfo(fs,fftsize); % frequency of each bin
%     fr_mask = fr_band(1) <= fr_full & fr_full <= fr_band(2); % mask for band
%     freqs = fr_full(fr_mask); % vector of just frequencies in band
%     ff_partial = ones(size(freqs)); % delta function
%     ff_partial = ff_partial + .5*randn(size(freqs)); % add noise
%     ff_partial = ff_partial .* exp(-1i*2*pi*freqs*tau); % delay
%     ff_full = fftsynth(ff_partial,freqs(1),fr_full,2); % fourier synthesis
%     yy = real(ifft(ff_full)); % convert to time domain
%     tt = linspace(0,(fftsize-1)/fs,fftsize); % time axis
%     clf; plot(tt,yy); % plot


if ~exist('dim','var'), dim=1; end

if dim==2, spect_partial = spect_partial.'; end

num_chans = size(spect_partial,2);
% find the first bin spect_partial will fit into
frmin_idx = find(freq_low(1)==freqs_full,1,'first');
assert(~isempty(frmin_idx),...
    'freq_low does not match any frequencies in freqs_full');
% zero-pad frequencies leading up to first bin
spect_full = [ zeros(frmin_idx-1,num_chans) ; spect_partial ];
% find the nyquist bin. if freqs_full is even, nyquist bin is negative, so
% use abs
[~,frnyq_idx] = max(abs(freqs_full));
assert(size(spect_full,1) <= frnyq_idx,...
    'supplied spect_partial has more bins than the positive frequencies in freqs_full');
% zero-pad up to the nyquist bin (the first nyquist bin in the case
% freqs_full is of odd length
spect_full = [ spect_full ; zeros(frnyq_idx-size(spect_full,1),num_chans) ];
% all that should remain are the negative frequency bins
num_negative_bins = length(freqs_full) - size(spect_full,1);
% bin 0 is DC, so start at idx 2. for real signals, the negative bins are
% mirrored over the middle (nyquist) bin and conjugated.
spect_full = [ spect_full ; conj(spect_full(num_negative_bins+1:-1:2,:)) ];

if dim==2, spect_full = spect_full.'; end
